home *** CD-ROM | disk | FTP | other *** search
/ Delphi Magazine Collection 2001 / Delphi Magazine Collection 20001 (2001).iso / DISKS / Issue33 / random / TstRndU3.pas < prev   
Encoding:
Pascal/Delphi Source File  |  1998-03-18  |  9.6 KB  |  300 lines
  1. {*********************************************************}
  2. {* TstRndU3                                              *}
  3. {* Copyright (c) Julian M Bucknall 1998                  *}
  4. {* All rights reserved.                                  *}
  5. {*********************************************************}
  6. {* Random number test program - tests                    *}
  7. {*********************************************************}
  8.  
  9. {Note: this unit is released as freeware. In other words, you are free
  10.        to use this unit in your own applications, however I retain all
  11.        copyright to the code. JMB}
  12.        
  13. unit TstRndU3;
  14.  
  15. interface
  16.  
  17. uses
  18.   TstRndU2;
  19.  
  20. {$I CHITABLE.INC}
  21.  
  22. procedure UniformityTest(RandGen     : TRandomGenerator;
  23.                      var ChiSquare   : double;
  24.                      var DegsFreedom : integer);
  25.  
  26. procedure GapTest(RandGen      : TRandomGenerator;
  27.                   Lower, Upper : double;
  28.               var ChiSquare    : double;
  29.               var DegsFreedom  : integer);
  30.  
  31. procedure PokerTest(RandGen     : TRandomGenerator;
  32.                 var ChiSquare   : double;
  33.                 var DegsFreedom : integer);
  34.  
  35. procedure CouponCollectorsTest(RandGen     : TRandomGenerator;
  36.                            var ChiSquare   : double;
  37.                            var DegsFreedom : integer);
  38.                 
  39. implementation
  40.  
  41. const
  42.   UniformityCount = 10000;
  43.   UniformityIntervals = 100;
  44.   GapBucketCount = 10;
  45.   GapsCount = 10000;
  46.   PokerCount = 50000;
  47.   CouponCount = 10000;
  48.  
  49. function IntPower(X  : double; N : integer) : double;
  50.   begin
  51.     case N of
  52.       0 : IntPower := 1.0;
  53.       1 : IntPower := X;
  54.       2 : IntPower := Sqr(X);
  55.     else
  56.       if Odd(N) then
  57.         IntPower := Sqr(IntPower(X, N div 2)) * X
  58.       else
  59.         IntPower := Sqr(IntPower(X, N div 2));
  60.     end;{case}
  61.   end;
  62.  
  63. function Stirling(N, M : integer) : double;
  64.   {-Calculate type-2 Stirling number, n & m are both assumed >= 0}
  65.   var
  66.     Work1, Work2 : double;
  67.   begin
  68.     if (N < M) then
  69.       Stirling := 0.0
  70.     else if (N = M) then
  71.       Stirling := 1.0
  72.     else
  73.       case m of
  74.         0 : Stirling := 0.0;
  75.         1 : Stirling := 1.0;
  76.         2 : Stirling := IntPower(2.0, n-1) - 1.0;
  77.       else
  78.         Work1 := Stirling(n-1, m);
  79.         Work2 := Stirling(n-1, m-1);
  80.         Stirling := (Work1 * m) + Work2;
  81.       end;{case}
  82.   end;
  83.  
  84. {=====================================================================
  85. The uniformity test.
  86. The random numbers are partitioned into a number of equally sized
  87. buckets between 0.0 and 1.0. On the average, each bucket should have
  88. the same number of random numbers; ie they should be evenly spread
  89. over the range [0.0, 0.1). Apply Chi-Squared test to the buckets.
  90. ======================================================================}
  91. procedure UniformityTest(RandGen     : TRandomGenerator;
  92.                      var ChiSquare   : double;
  93.                      var DegsFreedom : integer);
  94. var
  95.   BucketNumber,
  96.   i : integer;
  97.   Expected, ChiSqVal : double;
  98.   Bucket : array [0..pred(UniformityIntervals)] of integer;
  99. begin
  100.   {Fill buckets}
  101.   FillChar(Bucket, sizeof(Bucket), 0);
  102.   for i := 0 to pred(UniformityCount) do begin
  103.     BucketNumber := trunc(RandGen * UniformityIntervals);
  104.     inc(Bucket[BucketNumber]);
  105.   end;
  106.   {calc chi squared}
  107.   Expected := UniformityCount / UniformityIntervals;
  108.   ChiSqVal := 0.0;
  109.   for i := 0 to pred(UniformityIntervals) do
  110.     ChiSqVal := ChiSqVal + (Sqr(Expected - Bucket[i]) / Expected);
  111.   {return values}
  112.   ChiSquare := ChiSqVal;
  113.   DegsFreedom := pred(UniformityIntervals);
  114. end;
  115.  
  116. {=====================================================================
  117. The gap test.
  118. Each random number is tested to be in the range Lower..Upper. If it
  119. is a value of 1 is assigned, if not 0 is assigned. You'll get a stream
  120. of 0's and 1's. The lengths of the runs of 0's are then counted. These
  121. lengths are then bucketed, you'll get lengths of 0 upwards. These
  122. lengths are the 'gaps' between 1's. Apply Chi-Squared test to the
  123. buckets.
  124. ======================================================================}
  125. procedure GapTest(RandGen      : TRandomGenerator;
  126.                   Lower, Upper : double;
  127.               var ChiSquare    : double;
  128.               var DegsFreedom  : integer);
  129. var
  130.   NumGaps  : integer;
  131.   GapLen   : integer;
  132.   i        : integer;
  133.   p        : double;
  134.   Expected : double;
  135.   ChiSqVal : double;
  136.   R        : double;
  137.   Bucket   : array [0..pred(GapBucketCount)] of integer;
  138. begin
  139.   {calc gaps and fill buckets}
  140.   FillChar(Bucket, sizeof(Bucket), 0);
  141.   GapLen := 0;
  142.   NumGaps := 0;
  143.   while (NumGaps < GapsCount) do begin
  144.     R := RandGen;
  145.     if (Lower <= R) and (R < Upper) then begin
  146.       if (GapLen >= GapBucketCount) then
  147.         GapLen := pred(GapBucketCount);
  148.       inc(Bucket[GapLen]);
  149.       inc(NumGaps);
  150.       GapLen := 0;
  151.     end
  152.     else
  153.       if (GapLen < GapBucketCount) then
  154.         inc(GapLen);
  155.   end;
  156.   p := Upper - Lower;
  157.   ChiSqVal := 0.0;
  158.   {do all but the last bucket}
  159.   for i := 0 to GapBucketCount-2 do begin
  160.     Expected := p * IntPower(1-p, i) * NumGaps;
  161.     ChiSqVal := ChiSqVal + (Sqr(Expected - Bucket[i]) / Expected);
  162.   end;
  163.   {do the last bucket}
  164.   i := pred(GapBucketCount);
  165.   Expected := IntPower(1-p, i) * NumGaps;
  166.   ChiSqVal := ChiSqVal + (Sqr(Expected - Bucket[i]) / Expected);
  167.   {return values}
  168.   ChiSquare := ChiSqVal;
  169.   DegsFreedom := pred(GapBucketCount);
  170. end;
  171.  
  172.  
  173. {=====================================================================
  174. The poker test.
  175. The random numbers are grouped into 'hands' of 5, and the numbers are
  176. converted into a digit from 0..9. The number of different digits in
  177. each hand is then counted (1..5), and this result is bucketed. Because
  178. the probability of only one digit repeated 5 times is so low, it is
  179. grouped into the 2-different-digit category. Apply Chi-Squared test to
  180. the buckets.
  181. ======================================================================}
  182. procedure PokerTest(RandGen     : TRandomGenerator;
  183.                 var ChiSquare   : double;
  184.                 var DegsFreedom : integer);
  185. var
  186.   i, j,
  187.   BucketNumber,
  188.   NumFives : integer;
  189.   Accum, Divisor,
  190.   Expected, ChiSqVal : double;
  191.   Bucket : array [0..4] of integer;
  192.   Flag : array [0..9] of boolean;
  193.   p : array [0..4] of double;
  194. begin
  195.   {prepare}
  196.   FillChar(Bucket, sizeof(Bucket), 0);
  197.   NumFives := PokerCount div 5;
  198.   {calc probabilities for each bucket, algorithm from Knuth}
  199.   Accum := 1.0;
  200.   Divisor := IntPower(10.0, 5);
  201.   for i := 0 to 4 do begin
  202.     Accum := Accum * (10.0 - i);
  203.     p[i] := Accum * Stirling(5, succ(i)) / Divisor;
  204.   end;
  205.   {for each group of five random numbers, convert all five to a
  206.    number between 1 and 10, count the number of different digits}
  207.   for i := 1 to NumFives do begin
  208.     FillChar(Flag, sizeof(Flag), 0);
  209.     for j := 1 to 5 do begin
  210.       Flag[trunc(RandGen * 10.0)] := true;
  211.     end;
  212.     BucketNumber := -1;
  213.     for j := 0 to 9 do
  214.       if Flag[j] then inc(BucketNumber);
  215.     inc(Bucket[BucketNumber]);
  216.   end;
  217.  
  218.   {Accumulate the first bucket into the second, do calc separately -
  219.    it'll be the sum of the 'all the same' and 'two different digits'
  220.    buckets}
  221.   inc(Bucket[1], Bucket[0]);
  222.   Expected := (p[0]+p[1]) * NumFives;
  223.   ChiSqVal := Sqr(Expected - Bucket[1]) / Expected;
  224.   {write the other buckets}
  225.   for i := 2 to 4 do begin
  226.     Expected := p[i] * NumFives;
  227.     ChiSqVal := ChiSqVal + (Sqr(Expected - Bucket[i]) / Expected);
  228.   end;
  229.   {return values}
  230.   ChiSquare := ChiSqVal;
  231.   DegsFreedom := 3;
  232. end;
  233.  
  234.  
  235. {=====================================================================
  236. The coupon collectors test.
  237. The random numbers are read one by one, converted into a number from
  238. 0 to 4. The length of the sequence required to get a complete set of
  239. the digits 0..4 is counted, this will vary from 5 upwards. Once a full
  240. set is obtained, start over. Bucket the lengths of these sequences.
  241. Apply Chi-Squared test to the buckets.
  242. ======================================================================}
  243. procedure CouponCollectorsTest(RandGen     : TRandomGenerator;
  244.                            var ChiSquare   : double;
  245.                            var DegsFreedom : integer);
  246. var
  247.   NumSeqs, LenSeq, NumVals, NewVal,
  248.   i : integer;
  249.   Expected, ChiSqVal : double;
  250.   Bucket : array [5..20] of integer;
  251.   Occurs : array [0..4] of boolean;
  252.   p : array [5..20] of double;
  253. begin
  254.   {calc probabilities for each bucket, algorithm from Knuth}
  255.   p[20] := 1.0;
  256.   for i := 5 to 19 do begin
  257.     p[i] := (120.0 * Stirling(i-1, 4)) / IntPower(5.0, i);
  258.     p[20] := p[20] - p[i];
  259.   end;
  260.   {an alternative to calculate the last probability value:
  261.     p[last] := 1.0 - ((120.0 * Stirling(last-1, 5)) / IntPower(5.0, last-1)); }
  262.  
  263.   NumSeqs := 0;
  264.   FillChar(Bucket, sizeof(Bucket), 0);
  265.   while (NumSeqs < CouponCount) do begin
  266.     {keep getting coupons (ie random numbers) until we have collected
  267.      all five}
  268.     LenSeq := 0;
  269.     NumVals := 0;
  270.     FillChar(Occurs, sizeof(Occurs), 0);
  271.     repeat
  272.       inc(LenSeq);
  273.       NewVal := trunc(RandGen * 5);
  274.       if not Occurs[NewVal] then begin
  275.         Occurs[NewVal] := true;
  276.         inc(NumVals);
  277.       end;
  278.     until (NumVals = 5);
  279.     {update the relevant bucket depending on the number of coupoms we
  280.      had to collect}
  281.     if (LenSeq > 20) then
  282.       LenSeq := 20;
  283.     inc(Bucket[LenSeq]);
  284.     inc(NumSeqs);
  285.   end;
  286.  
  287.   {calculate chisquare value}
  288.   ChiSqVal := 0.0;
  289.   for i := 5 to 20 do begin
  290.     Expected := p[i] * NumSeqs;
  291.     ChiSqVal := ChiSqVal + (Sqr(Expected - Bucket[i]) / Expected);
  292.   end;
  293.  
  294.   {return values}
  295.   ChiSquare := ChiSqVal;
  296.   DegsFreedom := 15;
  297. end;
  298.  
  299. end.
  300.